library(foreign)
library(ggplot2)
library(grid)
library(gridExtra)
library(scales)

rm(list=ls())

setwd("Yourpath")

source("fine_grid.R") # ggplot layers
source("summarySE.R") # summary stats (for figure 5)


#################################################
#                                               #
#      PLOT TURNOUT PER CANTON  (Figure 5)      #
#                                               #
#################################################


# (i) Read in data

data = read.dta("data1.dta")
data <- data[data$canton!="AG" &data$canton!="AI" &data$canton!="AR"   &data$canton!="GL" &data$canton!="GR" &data$canton!="SG" &data$canton!="SH"  &data$canton!="TG" &data$canton!="TI" &data$canton!="ZH",]
X=read.csv("Y_treated_nonannual.txt", sep = "\t")[,1] 


# (ii) Compute turnout per referendumd day

dfc <- summarySE(data, measurevar="turnout_mean", groupvars=c("canton","dat_num"))



# (iii) Equally spaced

yearcutoffs <- c(1,11,21,42,59,74,106,127,130)

X_equalspace <- matrix(NA,yearcutoffs[9],1)

for (i in 2:length(yearcutoffs)){
  X_min <- min(X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])])
  X_max <- max(X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])])
  X_dist <- X_max-X_min 
  X_equalspace[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])]  <- (i-1)+ (X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])]-X_min)/X_dist 
}

X_equal_10years <- X_equalspace[is.element(X,yearcutoffs)]
X_equalspace[c(128:130)] <- (length(yearcutoffs)-1)+(X_equalspace[127]-X_equalspace[126])*c(1:3)

dfc$X_equalspace <- rep(X_equalspace,length(levels(factor(data$canton))))


x_breaks <- dfc$X_equalspace[is.element(dfc$dat_num,yearcutoffs) & dfc$canton=="VD"][1:8]
x_labels=seq(1900,1970,10)

dfc$canton_num <- as.numeric(as.factor(dfc$canton))
dfc$canton_num <- dfc$canton_num+1
dfc$canton_num[dfc$canton=="VD"] <- 1
dfc$canton_num[dfc$canton=="SZ"] <- dfc$canton_num[dfc$canton=="SZ"]-1
dfc$canton_num[dfc$canton=="SO"] <- dfc$canton_num[dfc$canton=="SO"]+1

dfc <- dfc[order(dfc$canton_num),] 

dfc$canton_num <- as.factor(dfc$canton_num)
levels_of_canton <- c(dfc$canton[dfc$X_equalspace==1& dfc$canton=="VD"],dfc$canton[dfc$X_equalspace==1& dfc$canton!="VD"])
levels_of_canton[1]<- "VD (treated)"
levels(dfc$canton_num) <- levels_of_canton

dfc$colour <- 0
dfc$colour[dfc$canton=="VD (treated)"] <- 1


# (iv) Plot

ggplot(dfc, aes(y=turnout_mean,x=X_equalspace,group(kkurz))) +
  geom_point(aes(shape=factor(colour))) + 
  geom_line() + 
  facet_wrap( ~canton_num, ncol = 4) +
  theme_bw_finegrid(base_size=26) +
  ylab("Turnout\n") + xlab("\n Year") + 
  scale_y_continuous(lim=c(0,1.02),labels=percent) +
  scale_x_continuous(breaks=x_breaks,labels=x_labels) +
  scale_shape_manual(values=c(16,1)) +
  theme(panel.margin = unit(1.6, "lines")) +
  theme(axis.text.x = element_text(size=14)) +
  theme(axis.text.y = element_text(size=16)) +
  guides(shape=FALSE) 
  

ggsave(file="Figure6.pdf",width=17,height=10 )



#################################################################
#                                                               #
#      PLOT RESULTS EXCL. PRE-TREATMENT OUTCOME  (Figure 6)     #
#                                                               #
#################################################################




Y_0 = read.csv("Y_synthetic_nonannual_cov_only.txt", sep = "\t")[,2]*100
Y_1 = read.csv("Y_treated_nonannual_cov_only.txt", sep = "\t")[,2]*100
X=read.csv("Y_treated_nonannual_cov_only.txt", sep = "\t")[,1] 


yearcutoffs <- c(1,11,21,42,59,74,106,127,130)

X_equalspace <- matrix(NA,length(Y_0),1)

for (i in 2:length(yearcutoffs)){
  X_min <- min(X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])])
  X_max <- max(X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])])
  X_dist <- X_max-X_min 
  X_equalspace[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])]  <- (i-1)+ (X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])]-X_min)/X_dist 
}

X_equal_10years <- X_equalspace[is.element(X,yearcutoffs)]
X_equalspace[c(128:130)] <- (length(yearcutoffs)-1)+(X_equalspace[127]-X_equalspace[126])*c(1:3)




gray.scale = .7
gray.scale2 = .7
dim.factor = .1
Legend = c("Vaud                  ","Synthetic Vaud")
Ylab = c("Turnout (%)")
Xlab = c("Year") 
Ylim = c(0,100) 
Xlim = c(min(X_equalspace),max(X_equalspace))
gap <- Y_1 - Y_0
plot(X_equalspace, Y_1, t = "l", col = "white", 
     lwd = 1,  ylab = Ylab, xlab = Xlab, ylim = Ylim, xlim=Xlim,
     yaxs = "i",xaxt="n") 
polygon(x=c(3.571429,3.571429,4.941176,4.941176),
        y=c(Ylim[1]+dim.factor,Ylim[2]-dim.factor,Ylim[2]-dim.factor,Ylim[1]+dim.factor),
        col=gray(gray.scale),border=gray(gray.scale),lty=1, bg=TRUE)
polygon(x=c( 5.466667, 5.466667,5.733333,5.733333),
        y=c(Ylim[1]+dim.factor,Ylim[2]-dim.factor,Ylim[2]-dim.factor,Ylim[1]+dim.factor),
        col=gray(gray.scale),border=gray(gray.scale),lty=1, bg=TRUE) 
lines(X_equalspace, Y_1, lwd = 1)
lines(X_equalspace, Y_0, col = gray(.55), lwd = 1)
points(X_equalspace, Y_1, pch=19, cex = 4/7 )
points(X_equalspace, Y_0, pch=19, col = gray(.55), cex = 4/7)

axis(1,las=1,labels=c(seq(1900,1970,10)),at=X_equal_10years[1:8])

legend(c("topright"), legend = Legend, lty = 1:1, col = c("black", 
                                                          gray(.55)), lwd = c(1.75, 1.75), cex = 6/7)

dev.copy2pdf(file="Figure7.pdf",width=8,height=6)


